Spatial Classes

Has slots for bbox which are the spatial extents, and a CRS object for the coordinate reference system.

getClass("Spatial")
## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2
getClass("CRS")
## Class "CRS" [package "sp"]
## 
## Slots:
##                 
## Name:   projargs
## Class: character
# Example instantiation
bb <- matrix(c(30, 85, 300, 90), ncol = 2, dimnames = list(NULL,c("min", "max")))
crs <- CRS("+proj=longlat")
Spatial(bb, crs)
## class       : Spatial 
## features    : 1 
## extent      : 30, 300, 85, 90  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

Spatial Points

Class SpatialPoints adds a coords slot to Spatial

getClass("SpatialPoints")
## Class "SpatialPoints" [package "sp"]
## 
## Slots:
##                                           
## Name:       coords        bbox proj4string
## Class:      matrix      matrix         CRS
## 
## Extends: "Spatial", "SpatialVector"
## 
## Known Subclasses: 
## Class "SpatialPointsDataFrame", directly
## Class "SpatialPixels", directly
## Class "SpatialPixelsDataFrame", by class "SpatialPixels", distance 2
# example instantiation
cran <- read.table("data/CRAN051001a.txt", header = TRUE)
cran_mat <- cbind(cran$long, cran$lat)

llcrs <- CRS("+proj=longlat +ellps=WGS84")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
## definition
cran_sp <- SpatialPoints(cran_mat, proj4string = llcrs)
summary(cran_sp)
## Object of class SpatialPoints
## Coordinates:
##                  min      max
## coords.x1 -122.95000 153.0333
## coords.x2  -37.81667  57.0500
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54
# accessor methods
coordinates(cran_sp)[which(cran$loc == "Brazil"),] # index coordinates
##      coords.x1 coords.x2
## [1,] -49.26667 -25.41667
## [2,] -42.86667 -20.75000
## [3,] -43.20000 -22.90000
## [4,] -47.63333 -22.71667
## [5,] -46.63333 -23.53333
bbox(cran_sp) # bounding box
##                  min      max
## coords.x1 -122.95000 153.0333
## coords.x2  -37.81667  57.0500
summary(cran_sp[which(cran$loc == "Brazil"), ]) # can index class directly, new bbox.
## Object of class SpatialPoints
## Coordinates:
##                 min       max
## coords.x1 -49.26667 -42.86667
## coords.x2 -25.41667 -20.75000
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 5

SpatialPointsDataFrame

Used for matching coordinates with other covariates, like name of coordinate. Needs to have matching row names. Objects are meant to operate like dataframes.

# matches by rownames
row.names(cran_mat) <- 1:nrow(cran)
# cran %>% row.names() # matches

cran_spdf <- SpatialPointsDataFrame(cran_mat, cran, proj4string = llcrs, match.ID = TRUE)
summary(cran_spdf) # meant to
## Object of class SpatialPointsDataFrame
## Coordinates:
##                  min      max
## coords.x1 -122.95000 153.0333
## coords.x2  -37.81667  57.0500
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54
## Data attributes:
##     place              north               east               loc           
##  Length:54          Length:54          Length:54          Length:54         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##       long                lat        
##  Min.   :-122.9500   Min.   :-37.82  
##  1st Qu.: -47.3833   1st Qu.: 34.52  
##  Median :   7.8500   Median : 42.73  
##  Mean   :  -0.6617   Mean   : 31.73  
##  3rd Qu.:  16.8333   3rd Qu.: 47.58  
##  Max.   : 153.0333   Max.   : 57.05
# alternate way, just attach coordinates and proj4string by assignment
cran_spdf_alt <- cran
# coordinates(cran_spdf_alt) <- c("long", "lat") # can select out columns if already in dataframe
coordinates(cran_spdf_alt) <- cran_mat # coordinates as matrix
proj4string(cran_spdf_alt) <- llcrs # crs
summary(cran_spdf_alt) # Turns it into a SpatialPointsDataFrame by assigning coordinate/proj
## Object of class SpatialPointsDataFrame
## Coordinates:
##                  min      max
## coords.x1 -122.95000 153.0333
## coords.x2  -37.81667  57.0500
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 54
## Data attributes:
##     place              north               east               loc           
##  Length:54          Length:54          Length:54          Length:54         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##       long                lat        
##  Min.   :-122.9500   Min.   :-37.82  
##  1st Qu.: -47.3833   1st Qu.: 34.52  
##  Median :   7.8500   Median : 42.73  
##  Mean   :  -0.6617   Mean   : 31.73  
##  3rd Qu.:  16.8333   3rd Qu.: 47.58  
##  Max.   : 153.0333   Max.   : 57.05

spacetime package

The spacetime package provides dataframes for certain types of spatial data.

  • STFDF
getClass("ST")
## Class "ST" [package "spacetime"]
## 
## Slots:
##                               
## Name:       sp    time endTime
## Class: Spatial     xts POSIXct
## 
## Known Subclasses: 
## Class "STF", directly
## Class "STS", directly
## Class "STI", directly
## Class "STT", directly
## Class "STFDF", by class "STF", distance 2
## Class "STSDF", by class "STS", distance 2
## Class "STIDF", by class "STI", distance 2
## Class "STTDF", by class "STT", distance 2
# Inheritance

Conditional AutoRegressive Model

Examples

Scotland Lip Cancer (Spatial CAR)

This section follows the tutorial from Geospatial Health Data: Modeling and Visualization with R-INLA and Shiny

data("scotland")

# getClass("SpatialPolygons") # has polygons, plotOrder, bbox, proj4string
# slotNames(map) # getting the value names of the object itself
map <- scotland$spatial.polygon
plot(map) # class SpatialPolygons

The map is in the projection OSGB 1936/British National Grid which has EPSG code 27700. The proj4 string of this projection can be seen in https://spatialreference.org/ref/epsg/27700/proj4/ or can be obtained with R as follows:

codes <- rgdal::make_EPSG()
codes[which(codes$code == "27700"), ]
##       code                           note
## 6288 27700 OSGB36 / British National Grid
##                                                                                                                  prj4
## 6288 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
##               prj_method
## 6288 Transverse Mercator
# assign 
proj4string(map) <- "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=km +no_defs"
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
## Warning in ReplProj4string(obj, CRS(value)): A new CRS was assigned to an object with an existing CRS:
## +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransform
# reproject into WGS84 for leaflet
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84 +no_defs"))

Data preparation

  • county: id of each country
  • Y: observed number of lip cancer cases in each county
  • E: expected number of lip cancer cases in each county
  • AFF: proportion of population in agriculture, fishing and forestry
  • SIR: standardized incidence ratio
# prepare the dataset and calculate SIR
d <- scotland$data[,c("county.names", "cases", "expected", "AFF")]
names(d) <- c("county", "Y", "E", "AFF")
d$SIR <- d$Y / d$E # observed over expected

# prepare SpatialPolygonsDataFrame, SpatialPolygons + Data basically
rownames(d) <- d$county # set rownames to county for matching later
# getClass("SpatialPolygonsDataFrame") # just adds "data" slot to SpatialPolygons
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE) # make Object for leaflet
l <- leaflet(map) %>% addTiles()

# domain gives possible values to be mapped, returns function
pal <- colorNumeric(palette = "YlOrRd", domain = map$SIR)

# Create the tooltip label
labels <- sprintf("<strong> %s </strong> <br/>
  Observed: %s <br/> Expected: %s <br/>
  AFF: %s <br/> SIR: %s",
  map$county, map$Y, round(map$E, 2),
  map$AFF, round(map$SIR, 2)
) %>%
  lapply(htmltools::HTML)

l %>%
  addPolygons(
    color = "grey", weight = 1, # the outlines
    fillColor = ~ pal(SIR), fillOpacity = 0.5, # how to fill the shapes
    highlightOptions = highlightOptions(weight = 4), # thicker outline on hover
    label = labels, # the tooltip label
    labelOptions = labelOptions( # some general options for the labels
      style = list(
        "font-weight" = "normal",
        padding = "3px 8px"
      ),
      textsize = "15px", direction = "auto"
    )
  ) %>%
  # legend in bottom right
  addLegend(
    pal = pal, values = ~SIR, opacity = 0.5,
    title = "SIR", position = "bottomright"
  )

For the modeling piece, we model with CAR,

\[ \begin{aligned} Y_i &\sim Poisson(E_i \theta_i), i = 1,\dots , n \\ \log(\theta_i) &= \beta_0 + \beta_1 \times AFF_i + u_i + v_i \end{aligned} \] - \(E_i\) is expected count - \(\theta_i\) is relative risk in area \(i\) - \(AFF_i\) is the fishing/farming covariate - \(u_i\) is spatial component with conditional guassian model, we assum \(u_i|\mathbf{u_{-i}} \sim N(\bar u_{\delta_i}, \frac{\sigma^2_u}{n_{\delta_i}})\) - \(v_i\) is unstructured spatial effect

# create neighborhood matrix from scottish data

nb <- poly2nb(map) # "nb" class

head(nb) # just a list of neighbors by node
## [[1]]
## [1]  5  9 19
## 
## [[2]]
## [1]  7 10
## 
## [[3]]
## [1] 12
## 
## [[4]]
## [1] 18 20 28
## 
## [[5]]
## [1]  1 12 19
## 
## [[6]]
## [1] 0
nb # printing shows some summary statistics of the dataset, as well as "disconnected regions" with no neighbors
## Neighbour list object:
## Number of regions: 56 
## Number of nonzero links: 234 
## Percentage nonzero weights: 7.461735 
## Average number of links: 4.178571 
## 3 regions with no links:
## orkney shetland western.isles
# nb2mat(nb, style="B",zero.policy=T) # convert to an adjacency matrix

Can calculate the Moran’s I, allows testing of spatial autocorrelation in the data, defined as

\[ \begin{equation} I = \frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}} \frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(x_i-\bar{x})(x_j-\bar{x})}{\sum_{i=1}^{n}(x_i - \bar{x})^2} \end{equation} \]

  • \(w_{ij}\) is weight of neighborhood incidence
  • \(x_i\) is covariate at area \(i\).
set.ZeroPolicyOption(TRUE) # setting globally, in argument it didn't work
## [1] FALSE
nb_W <- nb2listw(nb, style = "W")
nb_global <- Szero(nb_W) # get global sum of weights

# direct calculation of the statistic
moran(x =map$SIR, 
      listw = nb_W,
      n = length(map$SIR),
      S0 = nb_global)
## $I
## [1] 0.5251582
## 
## $K
## [1] 5.023736
moran.plot(map$SIR, nb_W) # upper right and lower left are positively correlated, slope is ~ approximately the moran statistic

moran.mc(map$SIR, nb_W, 100) # Lots of spatial correlation, randomization test.
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  map$SIR 
## weights: nb_W  
## number of simulations + 1: 101 
## 
## statistic = 0.49702, observed rank = 101, p-value = 0.009901
## alternative hypothesis: greater

Not quite sure how to interpret this plot,

# transfer the neighbors object to INLA
nb2INLA("output/map.adj", nb) # export to file, in INLA format
g <- inla.read.graph(filename = "output/map.adj") # read the .adj file for INLA
# indices for random effects u, v
map$idareau <- 1:nrow(map@data) # u
map$idareav <- 1:nrow(map@data) # v

# f() denotes random effect, "besag"
formula <- Y ~ AFF +
  f(idareau, model = "besag", graph = g, scale.model = TRUE) + # spatial dependence
  f(idareav, model = "iid") # unstructured noise

# fit the formual
res <- inla(formula, # formula from above
  family = "poisson", # log link
  data = map@data, # the data file
  E = E, # "Known component in the mean for the Poisson likelihoods defined as E exp(eta)" see docs
  control.predictor = list(compute = TRUE), # compute posteriors of predictions
  control.compute=list(return.marginals.predictor=TRUE) # marginals of each predictor as well, for exceedence probabilities
)

summary(res) # output
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
##    working.directory = working.directory, ", " silent = silent, inla.mode 
##    = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
##    .parent.frame)") 
## Time used:
##     Pre = 3.52, Running = 5.56, Post = 0.0319, Total = 9.12 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.305 0.120     -0.538   -0.305     -0.068   NA   0
## AFF          4.327 1.276      1.746    4.353      6.765   NA   0
## 
## Random effects:
##   Name     Model
##     idareau Besags ICAR model
##    idareav IID model
## 
## Model hyperparameters:
##                           mean       sd 0.025quant 0.5quant 0.975quant mode
## Precision for idareau     4.14     1.44       2.02     3.91       7.62   NA
## Precision for idareav 19639.46 19918.35    1375.78 13644.97   72579.48   NA
## 
## Marginal log-Likelihood:  -189.72 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

The results are pretty standard for interpretation, we can get marginals of the parameters (AFF) with smoothed (spline) marginals

# results has an output of marginals for each of the fixed effects, on some grid

# plot(res$marginals.fixed$AFF)
# lines(marginal$x, marginal$y)

marginal <- inla.smarginal(res$marginals.fixed$AFF) # smoothed marginals
marginal <- data.frame(marginal)
ggplot(marginal, aes(x = x, y = y)) + geom_line() +
  labs(x = expression(beta[1]), y = "Density") +
  geom_vline(xintercept = 0, col = "black") + theme_bw()

Each fitted value for each of the counties can be extracted from the summaries.

head(res$summary.fitted.values)
##                         mean        sd 0.025quant 0.5quant 0.975quant mode
## fitted.Predictor.01 4.965762 1.4504631   2.643878 4.790508   8.291851   NA
## fitted.Predictor.02 4.396906 0.6752405   3.172452 4.362819   5.817591   NA
## fitted.Predictor.03 3.620712 1.0169095   1.918926 3.523866   5.883690   NA
## fitted.Predictor.04 3.083339 0.8946611   1.627708 2.983424   5.113472   NA
## fitted.Predictor.05 3.329251 0.7500895   2.042165 3.267024   4.974026   NA
## fitted.Predictor.06 2.975147 0.9194491   1.490216 2.869251   5.070599   NA

final version of the leaflet map, with relative risk, and 95% CI for RR.

# Add to map df data
map$RR <- res$summary.fitted.values[, "mean"]
map$LL <- res$summary.fitted.values[, "0.025quant"]
map$UL <- res$summary.fitted.values[, "0.975quant"]

# Leaflet code 
pal <- colorNumeric(palette = "YlOrRd", domain = map$RR)

labels <- sprintf("<strong> %s </strong> <br/>
  Observed: %s <br/> Expected: %s <br/>
  AFF: %s <br/> SIR: %s <br/> RR: %s (%s, %s)",
  map$county, map$Y, round(map$E, 2),
  map$AFF, round(map$SIR, 2), round(map$RR, 2),
  round(map$LL, 2), round(map$UL, 2)
) %>% lapply(htmltools::HTML)

lRR <- leaflet(map) %>%
  addTiles() %>%
  addPolygons(
    color = "grey", weight = 1, fillColor = ~ pal(RR),
    fillOpacity = 0.5,
    highlightOptions = highlightOptions(weight = 4),
    label = labels,
    labelOptions = labelOptions(
      style =
        list(
          "font-weight" = "normal",
          padding = "3px 8px"
        ),
      textsize = "15px", direction = "auto"
    )
  ) %>%
  addLegend(
    pal = pal, values = ~RR, opacity = 0.5, title = "RR",
    position = "bottomright"
  )
lRR
# exceedence probabilities, cut offs for action
# P(\theta_i > c), probability that a specific area i has relative risk higher than some threshold.
sapply(res$marginals.fitted.values, FUN = function(marg){1 - inla.pmarginal(q = 2, marginal = marg)})
## fitted.Predictor.01 fitted.Predictor.02 fitted.Predictor.03 fitted.Predictor.04 
##        9.975709e-01        9.999996e-01        9.667358e-01        9.067867e-01 
## fitted.Predictor.05 fitted.Predictor.06 fitted.Predictor.07 fitted.Predictor.08 
##        9.793448e-01        8.680805e-01        9.780170e-01        4.596383e-01 
## fitted.Predictor.09 fitted.Predictor.10 fitted.Predictor.11 fitted.Predictor.12 
##        5.130059e-01        9.720600e-01        6.550928e-01        9.442041e-01 
## fitted.Predictor.13 fitted.Predictor.14 fitted.Predictor.15 fitted.Predictor.16 
##        7.405638e-01        5.351114e-01        4.747477e-01        6.258552e-01 
## fitted.Predictor.17 fitted.Predictor.18 fitted.Predictor.19 fitted.Predictor.20 
##        5.192967e-01        1.280056e-02        4.852896e-01        8.802111e-02 
## fitted.Predictor.21 fitted.Predictor.22 fitted.Predictor.23 fitted.Predictor.24 
##        3.173132e-02        2.622425e-02        1.326642e-02        6.028813e-05 
## fitted.Predictor.25 fitted.Predictor.26 fitted.Predictor.27 fitted.Predictor.28 
##        3.415062e-03        1.684688e-03        1.839832e-03        2.110785e-03 
## fitted.Predictor.29 fitted.Predictor.30 fitted.Predictor.31 fitted.Predictor.32 
##        1.553052e-03        5.385300e-05        2.241450e-04        2.026798e-01 
## fitted.Predictor.33 fitted.Predictor.34 fitted.Predictor.35 fitted.Predictor.36 
##        1.940328e-03        0.000000e+00        3.589546e-05        1.326963e-04 
## fitted.Predictor.37 fitted.Predictor.38 fitted.Predictor.39 fitted.Predictor.40 
##        4.311477e-05        0.000000e+00        5.081891e-03        4.865651e-06 
## fitted.Predictor.41 fitted.Predictor.42 fitted.Predictor.43 fitted.Predictor.44 
##        0.000000e+00        0.000000e+00        6.511447e-03        0.000000e+00 
## fitted.Predictor.45 fitted.Predictor.46 fitted.Predictor.47 fitted.Predictor.48 
##        0.000000e+00        0.000000e+00        0.000000e+00        0.000000e+00 
## fitted.Predictor.49 fitted.Predictor.50 fitted.Predictor.51 fitted.Predictor.52 
##        0.000000e+00        0.000000e+00        0.000000e+00        0.000000e+00 
## fitted.Predictor.53 fitted.Predictor.54 fitted.Predictor.55 fitted.Predictor.56 
##        0.000000e+00        0.000000e+00        2.773613e-03        5.010538e-04

Lung Cancer Iowa (spatio-temporal modeling)

Had to pull the data directly from the archived package, it seems to have been taken off CRAN because the contact author is no longer active. Link to Archive

dohio <- read.csv("data/dataohiocomplete.csv")
head(dohio)
##   county gender race year y     n  NAME
## 1      1      1    1 1968 6  8912 Adams
## 2      1      1    1 1969 5  9139 Adams
## 3      1      1    1 1970 8  9455 Adams
## 4      1      1    1 1971 5  9876 Adams
## 5      1      1    1 1972 8 10281 Adams
## 6      1      1    1 1973 5 10876 Adams
map <- as_Spatial(st_read("data/fe_2007_39_county", layer="fe_2007_39_county")) # read shape files, not quite, then convert to Spatial for sp methods
## Reading layer `fe_2007_39_county' from data source 
##   `/Users/mliou/Box Sync/gradwisc/reference/cheatsheets/mydocs/website/static/rmd/spatial/data/fe_2007_39_county' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 88 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.8203 ymin: 38.40342 xmax: -80.5182 ymax: 42.32713
## Geodetic CRS:  NAD83
plot(map)

# observed cases
d <- aggregate(
  x = dohio$y,
  by = list(county = dohio$NAME, year = dohio$year),
  FUN = sum
)
names(d) <- c("county", "year", "Y")
head(d)
##      county year  Y
## 1     Adams 1968  6
## 2     Allen 1968 32
## 3   Ashland 1968 15
## 4 Ashtabula 1968 27
## 5    Athens 1968 12
## 6  Auglaize 1968  7
# reorder correctly,
dohio <- dohio[order(
  dohio$county,
  dohio$year,
  dohio$gender,
  dohio$race
), ]

# check sorting
# dohio[1:20, ]

# 2 race x 2 gender
n.strata <- 4
E <- expected( # function needs sorted data
  population = dohio$n,
  cases = dohio$y,
  n.strata = n.strata
)
nyears <- length(unique(dohio$year)) # number of years
countiesE <- rep(unique(dohio$NAME), # county vector
                 each = nyears)

ncounties <- length(unique(dohio$NAME))
yearsE <- rep(unique(dohio$year), # rep by counties
              times = ncounties)

# putting it all together
dE <- data.frame(county = countiesE, 
                 year = yearsE,
                 E = E)
head(dE)
##   county year         E
## 1  Adams 1968  8.278660
## 2  Adams 1969  8.501767
## 3  Adams 1970  8.779313
## 4  Adams 1971  9.175276
## 5  Adams 1972  9.548736
## 6  Adams 1973 10.099777
# put it back into the main dataset
d <- merge(d, dE, by = c("county", "year"))
head(d)
##   county year  Y         E
## 1  Adams 1968  6  8.278660
## 2  Adams 1969  5  8.501767
## 3  Adams 1970  9  8.779313
## 4  Adams 1971  6  9.175276
## 5  Adams 1972 10  9.548736
## 6  Adams 1973  7 10.099777
# calc SIR
d$SIR <- d$Y / d$E
head(d)
##   county year  Y         E       SIR
## 1  Adams 1968  6  8.278660 0.7247549
## 2  Adams 1969  5  8.501767 0.5881130
## 3  Adams 1970  9  8.779313 1.0251372
## 4  Adams 1971  6  9.175276 0.6539313
## 5  Adams 1972 10  9.548736 1.0472591
## 6  Adams 1973  7 10.099777 0.6930846
# reshape into lengthwise file with rows as counties and columns as years
dw <- reshape(d,
  timevar = "year",
  idvar = "county",
  direction = "wide"
)

map <- merge(map, dw, by.x = "NAME", by.y = "county")

map_sf <- st_as_sf(map) # into sf for ggplot
map_sf <- gather(map_sf, year, SIR, paste0("SIR.", 1968:1988)) # longwise data now
map_sf$year <- as.integer(substring(map_sf$year, 5, 8)) # just get the year string

# make plot
ggplot(map_sf) + 
  geom_sf(aes(fill = SIR)) + # map the shapefiles 
  facet_wrap(~year, dir = "h", ncol = 7) +
  ggtitle("SIR") + theme_bw() +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank()
  ) +
  scale_fill_gradient2(
    midpoint = 1, low = "blue", mid = "white", high = "red"
  )

ggplot(d, aes(x = year, y = SIR, 
                   group = county, color = county)) +
  geom_line() + geom_point(size = 2) + theme_bw() +
  theme(legend.position = "none") +
  gghighlight(county == "Adams")
## Warning: Tried to calculate with group_by(), but the calculation failed.
## Falling back to ungrouped filter operation...
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## Warning: Tried to calculate with group_by(), but the calculation failed.
## Falling back to ungrouped filter operation...
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## label_key: county

Now we actually model the data with a spatio-temporal component

nb <- poly2nb(map) # convert to adjmat
nb2INLA("data/map_ohio.adj", nb)
g <- inla.read.graph(filename = "data/map_ohio.adj")

d$idarea <- as.numeric(as.factor(d$county))
d$idarea1 <- d$idarea # make copy of the idarea, becuase two terms with spatial random componenet, but in INLA, only one variable can be associated to the `f()` at a time.
d$idtime <- 1 + d$year - min(d$year) # index years/dates from 1

# for documentation of how this matches with model
# inla.doc(what = "bym")
formula <- Y ~ f(idarea, model = "bym", graph = g) + # spatial component
  f(idarea1, idtime, model = "iid") + # time x area interaction
  idtime # global trend across time

res <- inla(formula,
  family = "poisson", 
  data = d, 
  E = E,
  control.predictor = list(compute = TRUE)
)
summary(res)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
##    working.directory = working.directory, ", " silent = silent, inla.mode 
##    = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
##    .parent.frame)") 
## Time used:
##     Pre = 4.1, Running = 1.91, Post = 0.0435, Total = 6.05 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.530 0.024     -0.578   -0.530     -0.483   NA   0
## idtime       0.037 0.001      0.035    0.037      0.039   NA   0
## 
## Random effects:
##   Name     Model
##     idarea BYM model
##    idarea1 IID model
## 
## Model hyperparameters:
##                                              mean       sd 0.025quant 0.5quant
## Precision for idarea (iid component)        24.53     4.21      17.20    24.21
## Precision for idarea (spatial component)  1890.38  1772.12     121.01  1358.39
## Precision for idarea1                    49638.84 17544.85   24285.89 46629.35
##                                          0.975quant mode
## Precision for idarea (iid component)          33.75   NA
## Precision for idarea (spatial component)    6552.57   NA
## Precision for idarea1                      92463.65   NA
## 
## Marginal log-Likelihood:  -5940.95 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

It’s hard to interpret these

# plot with RR values
d$RR <- res$summary.fitted.values[, "mean"]
d$LL <- res$summary.fitted.values[, "0.025quant"]
d$UL <- res$summary.fitted.values[, "0.975quant"]

map_sf <- merge(
  map_sf, d,
  by.x = c("NAME", "year"),
  by.y = c("county", "year")
)

ggplot(map_sf) + geom_sf(aes(fill = RR)) +
  facet_wrap(~year, dir = "h", ncol = 7) +
  ggtitle("RR") + theme_bw() +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank()
  ) +
  scale_fill_gradient2(
    midpoint = 1, low = "blue", mid = "white", high = "red"
  )

Dynamic Spatio Temporal Statistical Models

This exposition follows Spatio-Temporal Statistics With R (Wikle, Zammit-Mangion, Cressie)

Examples

1D (IDE)

This example assumes the Itero-Differential Equation (IDE) is given, and we will show how to simulate from a given IDE. The ID we consider is

\[ \begin{aligned} Y_t(s) = \int_{D_s} m(s, x; \mathbf{\theta_p})Y_{t-1}(x)\,dx + \eta_t(s), \quad s,x\in D_s \end{aligned} \]

  • \(Y_t(s)\) is the spatio-temporal process at time \(t\),
  • \(m(s,x;\mathbf{\theta_p})\) - is called the transition kernel, and governs the dynamics of the process
  • \(\theta_p\) are the parameters, normally estimated, but we fix them for now.
  • \(\eta_t(\cdot)\) is a spatial process.
set.seed(1)

# create grid for simulation, (between 0,1)

ds <- 0.01 # step for spatial grid
s_grid <- seq(0, 1, by = ds) # spatial grid
N <- length(s_grid)
nT <- 201 # number of time points
t_grid <- 0:(nT-1) # time grid
st_grid <- expand.grid(s = s_grid, t = t_grid) # space-time grid
m <- function(x, s,thetap) {
 gamma <- thetap[1]             # amplitude
 l <- thetap[2]                 # length
 offset <- thetap[3]            # offset
 D <- outer(s + offset, x, "-") # displacements
 gamma * exp(-D^2/l)            # kernel eval.
}
# see how the parameterization differs for a few
thetap <- list()
thetap[[1]] <- c(40, 0.0002, 0) # narrow, centered.
thetap[[2]] <- c(5.75, 0.01, 0) # wider
thetap[[3]] <- c(8, 0.005, 0.1) # shifted right
thetap[[4]] <- c(8, 0.005, -0.1) # shifted left
m_x <- m(s=0.5, x = s_grid,
         thetap=thetap[[1]]) %>% 
  as.numeric()
df <- data.frame(x=s_grid, m = m_x)

# plot the data
ggplot(df) + 
  geom_line(aes(x,m)) +
  theme_bw() 

# The covariance of the spatial noise
# Simulate from Multivariate Normal, by cholesky decomp of the covariance
Sigma_eta <- 0.1 * exp( -abs(outer(s_grid, s_grid, '-') / 0.1))

L <- t(chol(Sigma_eta)) # chol() returns upper Cholesky factor
sim <- L %*% rnorm(nrow(Sigma_eta)) # simulate realization of eta

qplot(s_grid, sim, geom = "line") +
  theme_minimal() +
  labs(title = "realization of a random spatial process in 1d")

Y <- list()

# outer loop does 4 kernels
for(i in 1:4) {
  M <- m(s_grid, s_grid, thetap[[i]]) # create kernel function
  
  
  Y[[i]] <- data.frame(s = s_grid, # initialize grid and process
                       t = 0,
                       Y = 0)
  # for each timepoint
  for (j in t_grid[-1]) {
    prev_Y <- filter(Y[[i]], t == j - 1)$Y # get Y at t-1
    eta <- L %*% rnorm(N) # simulate eta
    
    new_Y <- (M %*% prev_Y * ds + eta) %>% as.numeric() # euler approximation, approximating integral with sum
    
    Y[[i]] <- rbind(Y[[i]],
                    data.frame(s = s_grid,
                               t = j,
                               Y = new_Y))
  }
}
hovmoller_plot <- function(idx) {
  ggplot(Y[[idx]]) + 
    geom_tile(aes(s,t,fill = Y)) +
    scale_x_continuous(expand = c(0,0)) +
    scale_y_reverse(expand = c(0,0)) + 
    fill_scale(name = "Y") +
    theme_bw() +
    theme(legend.position = "bottom")
}
p1 <- hovmoller_plot(1)
p2 <- hovmoller_plot(2)
p3 <- hovmoller_plot(3)
p4 <- hovmoller_plot(4)
plot_grid(p1, p2, p3, p4, nrow = 1)

Now we learn how to simulate data. The above graphs show the entire realization of the process model. In reality, we don’t get to observe all that data, instead, we just get little bits of the dynamics. So we create an “Incidence Matrix” \(\mathbf{H}_t\) to account for the data that we actually do observe. i.e., the above was the “process model”, now we implement the “data model”. We also assume we make some measurement errors, that are iid w/ mean zero.

\[ \begin{aligned} Z_t = H_tY_t + \varepsilon_t \end{aligned} \]

nobs <- 50 # number of observations at each timepoint
sobs <- sample(s_grid, nobs) # sample points on the grid
# should be w/ replacment right?

# Initialize incidence, nobs x N matrix
Ht <- matrix(0, nrow = nobs, ncol = N) # the incidence matrix
for(i in 1:nobs){
  idx <- which(sobs[i] == s_grid) # find elements and set to 1
  Ht[i, idx] <- 1
}
z_df <- data.frame() # init data frame to save noisy observations

# for each time point, we get the underlying process, take our observations and add measurement error
for (j in 0:(nT-1)) {
  Yt <- Y[[1]] %>% filter(t == j) %>% pull(Y) # get simulated process
  zt <- Ht %*% Yt + rnorm(nobs) # map 
  z_df <- rbind(z_df, data.frame(s = sobs, t = j, z = zt))
}
ggplot(z_df) + geom_point(aes(s, t, colour = z)) + 
  col_scale(name = "z") +
  scale_y_reverse(expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) + 
  theme_test() +
  labs(title = "Noisy Observation of Latent Process")

Spatio-Temporal Inference with the IDE Model

First example if simulating from a spatially invariant kernel, this is what k_spat_invariant=1 specificies. The movement here is very simple, it simply drifts away.

SIM1 <- simIDE(T = 10, nobs = 100, k_spat_invariant = 1)
str(SIM1, max.level=1)
## List of 6
##  $ s_df    :'data.frame':    16810 obs. of  4 variables:
##  $ z_df    :'data.frame':    1000 obs. of  4 variables:
##  $ z_STIDF :Formal class 'STIDF' [package "spacetime"] with 4 slots
##  $ g_truth :List of 9
##   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
##  $ g_obs   :List of 9
##   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
##  $ IDEmodel:List of 5
##   ..- attr(*, "class")= chr "IDE"
  • s_df - is the space df
  • z_df - is the observed noisy data
  • g_truth - a ggplot of the true process
  • g_obs - a ggplot of the observations
print(SIM1$g_truth)

print(SIM1$g_obs)

# returns object of class IDE and initial parameter estimates of \alpha_t
IDEmodel <- IDE(f=z~s1+s2,
                data = SIM1$z_STIDF,
                dt = as.difftime(1, units = "days"),
                grid_size = 41)

# Fits with deOptim, a derivative free optimzation method (same category as Nelder Mead and Simulated Annealing, but with different principles)
# Stands for differential evolution
# function is incredibly intenstive for calculation so the results can also be loaded below - 
# fit_results_sim1 <- fit.IDE(IDEmodel, parallelType = 1)
data("IDE_Sim1_results", package = "STRbook")
show_kernel(fit_results_sim1$IDEmodel) # show the kernel of the function
## Kernel is spatially invariant, plotting it centred on the origin.

# parameters of the model can be extracted from the object here
# real values are c(150, .002, -0.1, 0.1)
fit_results_sim1$IDEmodel$get("k") %>% unlist()
##          par1          par2          par3          par4 
## 152.836345912   0.001977115  -0.101601099   0.100368743
coef(fit_results_sim1$IDEmodel) # true values are c(.2, .2, .2)
## Intercept        s1        s2 
## 0.2073442 0.1966224 0.1907062

Now we extract the eigenvalues of the fitted evolution matrix to understand its long term behavior.

abs_ev <- eigen(fit_results_sim1$IDEmodel$get("M"))$values %>% 
  abs()

abs_ev %>% head() 
## [1] 0.4640858 0.4640858 0.4492270 0.4492270 0.4289470 0.4289470

The largest of these values is greater than 1, thus we conclude that the dynamics of this evolution are stable.

ST_grid_df <- predict(fit_results_sim1$IDEmodel) # can also take a new prediction grid, but by default it takes the simulation grid that we fit the data model for.
# Prediction map
gpred <- ggplot(ST_grid_df) + # Plot the predictions geom_tile(aes(s1, s2, fill=Ypred)) +
  geom_tile(aes(s1,s2, fill = Ypred)) +
  facet_wrap(~t) +
  fill_scale(name = "Ypred", limits = c(-0.1, 1.4)) + 
  coord_fixed(xlim=c(0, 1), ylim = c(0, 1)) +
  labs(title = "Predictions")

# Prediction standard error map
gpredse <- ggplot(ST_grid_df) + # Plot the prediction s.es
  geom_tile(aes(s1, s2, fill = Ypredse)) +
  facet_wrap(~t) +
  fill_scale(name = "Ypredse") +
  coord_fixed(xlim=c(0, 1), ylim = c(0, 1)) +
  labs(title = "Predictions SE")

plot_grid(gpred, gpredse)

### Now we run an example of spatially varying kernel

spatially varying kerneles are necessary in order to model forms of advection. specified with k_spat_invariant = 0

SIM2 <- simIDE(T = 15, nobs = 1000, k_spat_invariant = 0)
# the two ggplots of the truth and the simulated observations
print(SIM2$g_truth)

print(SIM2$g_obs)

appears to counter clockwise rotate, and then come to standstill in bottom of the domain. The varying advection that generated this field is visualized as follows…

show_kernel(SIM2$IDEmodel, scale = 0.2) # scale says change the arrow sizes.
## Kernel is spatially variant, plotting displacements

#  bisquare basis functions, same class as used by FRK.
mbasis_1 <- auto_basis(manifold = plane(), #fns, on the plane
                       data = SIM2$z_STIDF, # data
                       nres = 1, # 1 resolution
                       type = 'bisquare')
show_basis(mbasis_1)
## Note: show_basis assumes spherical distance functions when plotting

# now in our kernel, we can assume that amplitude (thetam1) and the scale (thetam2) are both constant
kernel_basis <- list(thetam1 = constant_basis(),
                     thetam2 = constant_basis(),
                     thetam3 = mbasis_1,
                     thetam4 = mbasis_1)

# specify the model
IDEmodel <- IDE(f = z ~ s1 + s2 + 1,
                data = SIM2$z_STIDF,
                dt = as.difftime(1, units = "days"),
                grid_size = 41,
                kernel_basis = kernel_basis)

# fit the IDE, computationally intensive so load from cache
# fit_results_sim2 <- fit.IDE(IDEmodel, parallelType = 1,
#                            itermax = 400)
data("IDE_Sim2_results", package = "STRbook")

show_kernel(fit_results_sim2$IDEmodel)
## Kernel is spatially variant, plotting displacements

ST_grid_df2 <- predict(fit_results_sim2$IDEmodel)
gpred <- ggplot(ST_grid_df2) + # Plot the predictions geom_tile(aes(s1, s2, fill=Ypred)) +
  geom_tile(aes(s1, s2, fill = Ypred)) +
  facet_wrap(~t) +
  fill_scale(name = "Ypred", limits = c(-0.1, 1.4)) + 
  coord_fixed(xlim=c(0, 1), ylim = c(0, 1)) +
  labs(title = "Predictions")

These are the predictions that were

Sydney Radar Data Set

Above we did a simulation example, but here we do a real data example

data("radar_STIDF", package = "STRbook")
  • hindcast how many time intervals in to past to estimate
  • forecase how many time intervals into the future to estimate.
IDEmodel <- IDE(f = z~1,
                data = radar_STIDF,
                dt=as.difftime(10, units = "mins"),
                grid_size = 41,
                forecast = 2,
                hindcast = 2)

# fitmodel
# fit_results_radar <- fit.IDE(IDEmodel, parallelType = 1)
data("IDE_Radar_results", package = "STRbook") # load cached results
show_kernel(fit_results_radar$IDEmodel) 
## Kernel is spatially invariant, plotting it centred on the origin.

shifted off center, so suggestive of predominatly north east direction.

shift_pars <- (fit_results_radar$IDEmodel$get("k") %>% unlist())[3:4]
print(shift_pars)
##      par3      par4 
## -5.488385 -1.860784

We can make a calculation here, \(\sqrt{5.5^2 + 1.9^2} = 5.82\) in 10 minutes, 34.91 km per hour.

abs_ev <- eigen(fit_results_radar$IDEmodel$get("M"))$values %>% abs()
summary(abs_ev)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.008377 0.581449 0.671552 0.615239 0.719497 0.786783

maximum eigenvalue less than 1 means that its stable, but more persistent than the simulation example.

ST_grid_df <- predict(fit_results_radar$IDEmodel)
radar_df$time <- format(radar_df$t, "%H:%M")
ST_grid_df$time <- format(ST_grid_df$t, "%H:%M")
## Add time records with missing data
radar_df <- rbind.fill(radar_df, data.frame(time = c("08:05", "08:15",
                                           "10:25", "10:35")))

## Plot of data, with color scale capped to (-20, 60)
gobs <- ggplot(radar_df) +
  geom_tile(aes(s1, s2, fill = pmin(pmax(z, -20), 60))) + 
  fill_scale(limits = c(-20, 60), name = "Z") + 
  facet_wrap(~time) + 
  coord_fixed() + 
  theme_bw() +
  labs(title = "Observed Sydney Data")

## Plot of predictions with color scale forced to (-20, 60)
gpred <- ggplot(ST_grid_df) +
  geom_tile(aes(s1, s2, fill = Ypred)) + facet_wrap(~time) + coord_fixed() + theme_bw() + fill_scale(limits = c(-20, 60), name = "Ypred") +
  labs(title = "Predicted Data, with 2 time periods before and after")

plot_grid(gobs, gpred)
## Warning: Removed 4 rows containing missing values (geom_tile).

Inference with Unknown Evolution Operator

We use EOFs are ideal basis functions for exploring the observed signal.

We fit two models, a traditional one in time-series framework with vector autoregression and use method of moments.

# trucate at April 1997 in order to forcast 6 months ahead

data("SSTlandmask", package = "STRbook")
data("SSTlonlat", package = "STRbook")
data("SSTdata", package = "STRbook")
delete_rows <- which(SSTlandmask == 1) # remove land values
SST_Oct97 <- SSTdata[-delete_rows, 334] # save Oct 1997 SSTs  
SSTdata <- SSTdata[-delete_rows, 1:328] # until April 1997 , 399 times points normally, but truncate
SSTlonlat$mask <- SSTlandmask # assign mask to df
# 2520 x 328 dataframe
SSTdata %>% head()
##           V1        V2        V3        V4          V5          V6         V7
## 16 0.2296104 0.2515602 0.4510937 0.2989082 -0.28913879 -0.20289230 -0.1392193
## 17 0.4862499 0.4507828 0.5332813 0.5499210 -0.22593880 -0.05968857  0.2024994
## 18 0.7855454 0.6135120 0.5860958 0.6160946 -0.11890602 -0.01961136  0.3596077
## 19 0.9411717 0.7603111 0.7017193 0.5956268 -0.02531242 -0.02273369  0.4096889
## 20 1.0239792 0.7840633 0.7537537 0.6158600  0.02359390 -0.06570435  0.3028908
## 21 1.0313301 0.8165626 0.8502331 0.6375790  0.09546661 -0.03921890  0.2607803
##            V8          V9         V10         V11       V12       V13
## 16 -0.3317967 -0.08953095 0.000623703  0.16234589 0.3113289 0.6296101
## 17  0.0696888  0.15523529 0.237970350  0.11281204 0.5346890 0.7362499
## 18  0.2849217  0.20914078 0.395938870 -0.06546974 0.5482826 0.7855454
## 19  0.3793736  0.20929337 0.427658080 -0.05945587 0.5828152 0.8411713
## 20  0.4363270  0.26359367 0.372110370 -0.05234337 0.5239067 0.9239807
## 21  0.4079704  0.33273506 0.348985670 -0.01726532 0.4121094 0.9813290
##           V14         V15        V16         V17       V18        V19
## 16 -0.1484413 -0.29890633 -0.4010925 -0.18914032 0.1471081 -0.1392174
## 17  0.2507839 -0.16671944 -0.1000786  0.02406311 0.5903110  0.3024979
## 18  0.7135143 -0.06390381  0.1160946  0.18109322 0.9303894  0.5096092
## 19  1.0103111  0.05171967  0.2456264  0.32468605 1.0772667  0.5596886
## 20  1.0840626  0.10375214  0.3158608  0.47359467 1.0842972  0.5528908
## 21  1.0665646  0.20023537  0.3875771  0.69546700 1.1107807  0.5607834
##            V20        V21        V22         V23        V24        V25
## 16 -0.08179855 -0.2395306 -0.1993771 -0.58765602 -0.4886703 -0.3203907
## 17  0.06968880  0.1552353  0.1879711 -0.23718643 -0.3653126 -0.2137508
## 18  0.18492317  0.4591408  0.4959373 -0.01547051 -0.3517170 -0.3144531
## 19  0.22937393  0.6592922  0.6776581  0.14054489 -0.3171864 -0.5088291
## 20  0.28632736  0.7135944  0.6721096  0.14765739 -0.2760944 -0.6260185
## 21  0.30796814  0.6827354  0.6489849  0.08273506 -0.3378906 -0.6186695
##           V26        V27         V28         V29        V30         V31
## 16 -0.4984398 -0.5489063 -0.10109329 -0.53913879 -0.3528919 -0.08921814
## 17 -0.3492184 -0.6167183  0.14992332 -0.22593880 -0.2096882 -0.14750099
## 18 -0.3864861 -0.6639042  0.21609306  0.08109283 -0.1196098 -0.24039078
## 19 -0.4396877 -0.4982815  0.14562798  0.12468910 -0.1227322 -0.24031067
## 20 -0.6159363 -0.2462463  0.06586075  0.12359619 -0.2157040 -0.09710693
## 21 -0.6834374  0.1502342  0.08757782  0.09546661 -0.2392197  0.11078262
##           V32         V33       V34       V35         V36       V37       V38
## 16 0.06820297 -0.33953094 0.2006245 0.7123470 -0.03866959 0.3796101 0.7515602
## 17 0.21968651 -0.14476395 0.5379696 0.7128124 -0.31531143 0.5362492 0.9007816
## 18 0.13492393 -0.09085846 0.6959381 0.6345291 -0.50171661 0.6355476 0.9135132
## 19 0.07937241 -0.04070663 0.6276588 0.5905438 -0.41718483 0.6911717 0.8103104
## 20 0.03632736  0.06359291 0.4221096 0.5476570 -0.17609215 0.7239800 0.6840630
## 21 0.05797005  0.33273506 0.2989864 0.6327343  0.06211090 0.6813297 0.5665626
##          V39        V40         V41       V42       V43       V44       V45
## 16 0.7010937 -0.4510918 -0.23913956 0.7971077 0.7107811 0.7182026 0.9604683
## 17 0.7832794 -0.2000790 -0.02593803 0.8903122 0.8024979 0.7696876 1.0552330
## 18 0.7860966  0.1160946  0.23109245 0.8303890 0.7096100 0.6849232 0.9091415
## 19 0.7517204  0.3456287  0.52468681 0.6272678 0.6096878 0.6793728 0.7592945
## 20 0.8037529  0.4158592  0.72359467 0.5342960 0.5528908 0.6863270 0.6135922
## 21 0.8002338  0.3875771  0.79546738 0.5607796 0.4607811 0.6579704 0.5327339
##          V46       V47       V48       V49       V50        V51          V52
## 16 0.9506245 0.6123447 0.7113304 0.3296108 0.1515617 -0.5489063 -0.001091003
## 17 1.1879711 0.8128128 0.6846886 0.4362488 0.1507816 -0.5167198  0.049921036
## 18 1.2459373 0.9845295 0.5982838 0.4855461 0.1635151 -0.3139057  0.066095352
## 19 1.0776577 0.9905453 0.6328144 0.4911728 0.1603108 -0.1482792  0.145627980
## 20 1.0221100 0.9976559 0.7239056 0.5239792 0.1840630 -0.1462479  0.265859600
## 21 0.9489842 1.0327358 0.8121109 0.5813293 0.2665634 -0.1497650  0.387578960
##            V53          V54          V55         V56         V57        V58
## 16 -0.13913918 -0.102891920 -0.089218140 -0.28179550 -0.03952980 -0.3493748
## 17 -0.07593918 -0.009689331  0.002498627 -0.23031044  0.05523491 -0.2120304
## 18 -0.06890678 -0.169610980  0.009609222 -0.21507835 -0.04085922 -0.2540627
## 19 -0.12531281 -0.272733690  0.009689331 -0.12062645 -0.09070587 -0.2223415
## 20 -0.17640495 -0.165702820  0.102891920 -0.11367226 -0.13640785 -0.2778912
## 21 -0.05453491  0.060779572  0.310781480 -0.09202957 -0.21726608 -0.2510128
##           V59        V60         V61        V62         V63         V64
## 16 -0.8376541 -1.2386703 -0.77038956 -0.9984398 -0.14890671 -0.20109177
## 17 -0.7871876 -0.9153118 -0.56375122 -0.8492184  0.08328247  0.04992104
## 18 -0.7654705 -0.6017170 -0.41445351 -0.5864868  0.18609428  0.21609497
## 19 -0.6094551 -0.2671871 -0.30882835 -0.4896889  0.25172043  0.09562874
## 20 -0.3523426  0.0239048 -0.22602081 -0.3659363  0.25375366  0.06586075
## 21 -0.1672649  0.3121109 -0.01867104 -0.3334370  0.20023537  0.13757706
##          V65         V66       V67       V68       V69         V70          V71
## 16 0.1108608 -0.50289345 0.2607803 0.4182014 0.2604694 -0.34937477  0.212347030
## 17 0.3240623 -0.25968933 0.5024986 0.5196876 0.5052338 -0.21203041  0.212812420
## 18 0.3310947 -0.06961060 0.4596081 0.5849247 0.6091423 -0.05406189  0.134529110
## 19 0.2746868 -0.02273369 0.3096886 0.6793728 0.6092930  0.02765655  0.040544510
## 20 0.2235947  0.03429604 0.3028908 0.6863270 0.5635929  0.07210922 -0.002342224
## 21 0.2454681  0.16077995 0.4107819 0.6079712 0.4327354  0.09898567  0.082735062
##            V72         V73        V74        V75        V76        V77
## 16 -0.03866959 -0.17038918 -0.3984394 -0.1989078 -0.1510925 -0.2891388
## 17  0.18468857 -0.06375122 -0.6492176 -0.1167183 -0.3500767 -0.3259392
## 18  0.24828339 -0.01445389 -0.8364849 -0.1139050 -0.6339054 -0.4189053
## 19  0.18281364 -0.05882835 -0.9896908 -0.2482796 -0.8043728 -0.5253124
## 20  0.12390709 -0.12601852 -1.0659370 -0.2962456 -0.7841396 -0.6264038
## 21  0.21210861  0.03132820 -0.9334354 -0.3997669 -0.6624222 -0.6045341
##           V78         V79        V80        V81         V82        V83
## 16  0.2971077 -0.23921967 -0.4817963 -0.9395313  0.10062408 -0.1876545
## 17  0.1403122 -0.04750061 -0.4303112 -0.9447651 -0.01202965 -0.1871872
## 18 -0.0696106  0.05960846 -0.5650768 -0.9908581 -0.15406227 -0.2654705
## 19 -0.1727333  0.10968781 -0.6706276 -0.9907055 -0.27234268 -0.4594555
## 20 -0.2157040  0.00289154 -0.5636711 -0.9864063 -0.22789192 -0.6023426
## 21 -0.2392197 -0.33921814 -0.3420296 -0.9172649 -0.15101433 -0.7672653
##           V84         V85           V86        V87         V88         V89
## 16 -0.2386704  0.07961082  0.0015602112 -0.1989078  0.19890785  0.01086044
## 17 -0.1653118 -0.06374931  0.0007820129 -0.2167206  0.09992027 -0.07593918
## 18 -0.1017170 -0.26445389 -0.0864849090 -0.3139057 -0.13390541 -0.21890831
## 19 -0.1671848 -0.60882950 -0.2396888700 -0.4982815 -0.30437279 -0.32531166
## 20 -0.2260933 -0.97601891 -0.5159378100 -0.6962471 -0.43413925 -0.47640419
## 21 -0.3878899 -1.11866950 -0.6334362000 -0.8997669 -0.46242142 -0.55453300
##             V90        V91         V92         V93        V94        V95
## 16 -0.252893450 -0.6392193 -0.13179779  0.31046867 -0.1993771 -0.6376553
## 17 -0.009689331 -0.5475006  0.06968880  0.15523529 -0.1620293 -0.2371864
## 18  0.180389400 -0.6403923  0.08492279  0.05914116 -0.1540623  0.1345291
## 19  0.227266310 -0.6903114  0.07937241 -0.04070473 -0.3223419  0.3405437
## 20  0.134296420 -0.7971077  0.03632736 -0.13640785 -0.4778900  0.3476582
## 21  0.060779572 -0.9392185 -0.04203033 -0.31726456 -0.5010147  0.1327343
##           V96         V97        V98          V99         V100        V101
## 16 -0.4386711  0.12961006 0.35155869  0.001092911 -0.001092911  0.01086044
## 17 -0.4153118 -0.01375008 0.30078316  0.033281326 -0.100078580 -0.07593918
## 18 -0.3517170 -0.21445274 0.11351395 -0.063903809  0.016094208 -0.21890831
## 19 -0.3671856 -0.45882797 0.06031227 -0.298280720  0.045627594 -0.32531166
## 20 -0.4260941 -0.57601929 0.03406334 -0.496246340 -0.034141541 -0.37640381
## 21 -0.4878883 -0.76867104 0.01656151 -0.699766160 -0.112421040 -0.60453415
##          V102       V103       V104       V105       V106       V107
## 16 -0.1028919 -0.4892197 -0.3317967 -0.3895321 -0.4493771 -0.4876537
## 17 -0.1096878 -0.2475014 -0.2303123 -0.1447659 -0.6620293 -0.3871880
## 18 -0.2196102 -0.2403908 -0.5150776 -0.2408600 -0.8540630 -0.5154686
## 19 -0.3727322 -0.3903122 -0.7706261 -0.4407063 -1.0223427 -0.6594562
## 20 -0.5157051 -0.7471085 -1.0136719 -0.5364075 -1.1278896 -0.7523422
## 21 -0.5892181 -1.1392174 -1.1420307 -0.5672646 -1.2010136 -0.7672672
##           V108        V109        V110        V111       V112       V113
## 16 -0.08866882 -0.07039070 -0.29843903 -0.24890709 -0.2510929 0.06085968
## 17  0.08468819  0.03624916 -0.14921761 -0.11671829 -0.4000778 0.22406006
## 18  0.14828491 -0.01445389 -0.03648567 -0.01390266 -0.4839039 0.23109245
## 19 -0.06718636 -0.10882950  0.01031113  0.05171967 -0.5543728 0.17468834
## 20 -0.37609291 -0.12601852 -0.01593780  0.05375290 -0.4841404 0.17359543
## 21 -0.53788948 -0.06867027 -0.13343620  0.05023384 -0.4124222 0.14546585
##            V114      V115        V116       V117        V118      V119
## 16 -0.002893448 0.3607807 -0.03179741 -0.5395298 -0.09937668 0.2623463
## 17 -0.009689331 0.4024983  0.06968880 -0.3447647 -0.06202888 0.5128136
## 18 -0.219610210 0.2596092  0.08492279 -0.2408600 -0.25406265 0.6845303
## 19 -0.272733690 0.2096882  0.02937317 -0.1907062 -0.42234230 0.7405453
## 20 -0.265703200 0.1528912 -0.06367111 -0.1364078 -0.47789001 0.5976582
## 21 -0.139219280 0.1107826 -0.09202957 -0.1172638 -0.40101433 0.4827347
##         V120        V121        V122        V123        V124        V125
## 16 0.7113304  0.27961159 0.001560211  0.25109291  0.29890823 -0.03913879
## 17 1.0346889  0.53624916 0.100782390  0.13328171  0.29992104  0.02406311
## 18 1.0982838  0.58554840 0.163515090 -0.01390266  0.16609383 -0.11890602
## 19 1.0828152  0.54117203 0.160310750  0.05171967 -0.00437355 -0.42531204
## 20 1.0239048  0.22397995 0.134063720  0.10375214 -0.28413963 -0.57640457
## 21 0.7621098 -0.01867104 0.166563030  0.15023422 -0.51242065 -0.65453339
##           V126       V127      V128      V129      V130       V131        V132
## 16 -0.15289116 -0.5392189 0.3682022 0.5604687 0.5506229 0.46234512  0.36132812
## 17 -0.05968857 -0.6475029 0.4196873 0.6552353 0.9379711 0.61281204  0.18468857
## 18  0.03038979 -0.6403923 0.5849228 0.6591415 1.0959358 0.68452835  0.19828415
## 19  0.12726784 -0.4403114 0.5793724 0.6592941 1.1776581 0.54054451  0.03281403
## 20  0.18429565 -0.0471077 0.3363266 0.5635929 1.0721092 0.29765892 -0.27609444
## 21  0.16077995  0.3607826 0.1079712 0.2827339 0.8989849 0.03273583 -0.18788910
##          V133      V134      V135        V136        V137       V138       V139
## 16 -0.4703884 0.3015595 0.2010937  0.34890747  0.21086121 -0.1528931  0.2607803
## 17 -0.4137516 0.5007820 0.3832817  0.34992027  0.12406158 -0.1596889 -0.1475029
## 18 -0.2644539 0.6135139 0.4860954  0.11609459  0.03109169 -0.2696114 -0.6403923
## 19 -0.1588287 0.7603092 0.5517197 -0.00437355  0.02468681 -0.3227329 -0.7403107
## 20 -0.1260185 0.9340630 0.6037521  0.01585960 -0.07640457 -0.2657032 -0.4471092
## 21 -0.2186699 1.0165634 0.6002350  0.23757935  0.04546738 -0.1892204 -0.0392189
##          V140        V141        V142        V143      V144      V145
## 16 -0.7817974  0.16046906 -0.49937630  0.03609657 0.7950783 0.8483601
## 17 -0.9803104 -0.04476547 -0.41203117 -0.07968712 0.7746887 0.7725010
## 18 -1.0650768 -0.24085999 -0.35406303 -0.24671745 0.8395328 0.9817963
## 19 -1.0206261 -0.34070778 -0.22234154 -0.53445625 0.7090645 1.1211700
## 20 -0.7636719 -0.43640709 -0.07789040 -0.73484230 0.4301567 1.1452293
## 21 -0.5920296 -0.26726532  0.09898567 -0.64476585 0.1733589 1.0863323
##           V146        V147        V148      V149        V150        V151
## 16  1.06906130  0.77859306 0.637657170 0.9833584  0.33085823 -0.02671814
## 17  0.60828018  0.42952919 0.227422710 0.4353123 -0.08218956 -0.30250359
## 18  0.39851379  0.40234566 0.208593370 0.3185940 -0.21085930 -0.16664124
## 19  0.17281151  0.31171989 0.316875460 0.3934364 -0.16648483  0.02218819
## 20  0.01656342  0.10500336 0.229610440 0.3823452 -0.15695190  0.14789200
## 21 -0.05843735 -0.05351639 0.005079269 0.2579651 -0.16296959  0.20453262
##           V152        V153       V154       V155        V156       V157
## 16  0.12570190 -0.43577957 -0.8331261 -0.7339058 -0.21742058 -0.7816391
## 17 -0.16781044 -0.68476486 -0.9507809 -0.9159355 -0.37906075 -0.8950024
## 18 -0.09757805 -0.56586075 -0.8340626 -0.8817196 -0.29546738 -0.8582020
## 19  0.09687424 -0.30320549 -0.5198422 -0.6694565 -0.03718758 -0.5488319
## 20  0.26757812 -0.10640717 -0.3866406 -0.4398441  0.32140541 -0.2497692
## 21  0.23421860 -0.08226585 -0.4585152 -0.3460140  0.36960983 -0.2899208
##           V158      V159       V160      V161      V162        V163        V164
## 16 -0.30718994 0.1573410 0.05515671 0.3821087 0.6333599  0.91703033  0.61695290
## 17 -0.38046646 0.1532784 0.15867233 0.3565579 0.4078121  0.51374817  0.09093666
## 18 -0.33023643 0.2698460 0.40234566 0.3273430 0.3666382  0.45085907  0.06742287
## 19 -0.20718956 0.2617188 0.58687592 0.2296867 0.2435169  0.09093857 -0.08062744
## 20 -0.00718689 0.2062511 0.67585945 0.3248463 0.2292957 -0.14710999 -0.21367264
## 21 -0.02343750 0.1427345 0.61632729 0.3792152 0.1582813 -0.12171745 -0.06077957
##         V165        V166       V167        V168       V169       V170
## 16 0.7117176  1.38562390  1.0110931 -0.05742073 -0.5178909 -0.1209412
## 17 0.2564831  0.51922035  0.5528126 -0.60531235 -1.1012516 -0.3117180
## 18 0.3141422  0.10843658  0.1707783 -1.01671600 -1.3957024 -0.6764889
## 19 0.2755413 -0.09609222 -0.1082058 -1.14843750 -1.2488289 -0.8209400
## 20 0.1223450 -0.07289123 -0.1173420 -1.05734440 -1.0360203 -0.7784386
## 21 0.1727352  0.06648445  0.1189823 -0.86164093 -1.0849228 -0.7984371
##          V171       V172       V173        V174       V175       V176
## 16 -0.1339092 -0.4660950 -0.1666393 -0.08039475 0.12578011 0.21570206
## 17 -0.2142200 -0.5988293 -0.2109375 -0.24968910 0.06374741 0.04093742
## 18 -0.4201546 -0.6976547 -0.2576561 -0.20336151 0.23960876 0.02742195
## 19 -0.5032806 -0.5843735 -0.2228146 -0.16773415 0.29718781 0.07187271
## 20 -0.5249977 -0.5328903 -0.1551533 -0.05320358 0.33039093 0.16257668
## 21 -0.4910164 -0.5149231 -0.1007843 -0.12796974 0.20953178 0.11421967
##           V177       V178        V179       V180       V181       V182
## 16 -0.22078323 -0.5731258  0.12359619 -0.2949219  0.1958599 -0.1096897
## 17 -0.55101585 -0.8070278 -0.11468697 -0.3703117 -0.2112522 -0.5754681
## 18 -0.51210785 -0.7140636  0.04577828 -0.3154678 -0.5044517 -0.7539864
## 19 -0.30445671 -0.4910946  0.26179314 -0.2834358 -0.6175804 -0.7184391
## 20 -0.08765793 -0.3266411  0.32390594 -0.2535934 -0.4972687 -0.4896889
## 21 -0.11726379 -0.3822651  0.28148270 -0.1966400 -0.1874199 -0.3184357
##          V183       V184        V185       V186       V187        V188
## 16 -0.2151585  0.1864052  0.12085915 -0.2991428 -0.2454681 -0.07679748
## 17 -0.5679703 -0.1863270 -0.26968765 -0.8096905 -0.7050018 -0.30031395
## 18 -0.7351551 -0.3064079 -0.31140900 -0.7671108 -0.4891434 -0.26632690
## 19 -0.7907810 -0.2668724 -0.09656334 -0.3864842 -0.1265602 -0.24062729
## 20 -0.5812473 -0.1128921  0.01484489 -0.2807026 -0.1658592 -0.20242310
## 21 -0.4510174 -0.2074223 -0.18203354 -0.5467205 -0.4217186 -0.13328171
##          V189        V190       V191        V192      V193      V194       V195
## 16  0.2692184  0.36312485 -0.1226540  0.25008011 0.4158592 0.1753082 -0.4701595
## 17 -0.1335163 -0.11702919 -0.5834370  0.04843903 0.3225002 0.2020321 -0.2867184
## 18 -0.2958603 -0.41156387 -0.7179718 -0.14296722 0.2180481 0.3510132 -0.1589031
## 19 -0.4294567 -0.42484093 -0.5457058 -0.15093803 0.2024193 0.4028091 -0.2270317
## 20 -0.3414059 -0.17914009 -0.2060928 -0.10359383 0.2714806 0.3728123 -0.3724976
## 21 -0.1660156 -0.06351471 -0.1697655 -0.11788940 0.3925800 0.4190617 -0.2622662
##           V196      V197      V198       V199          V200        V201
## 16  0.07890701 0.6096096 1.0221043 0.47328186  0.4857006100  0.58547020
## 17  0.26742172 0.7765636 0.7503109 0.12499809 -0.0790615080  0.18398285
## 18  0.19734573 0.6723423 0.4016399 0.03335762 -0.1675758400  0.01664162
## 19 -0.10312271 0.3846874 0.1810169 0.13718987  0.0006217956 -0.07320786
## 20 -0.36289024 0.2610950 0.2867947 0.21664047  0.0688266750 -0.10890579
## 21 -0.30367279 0.4079666 0.4595299 0.26703262  0.0092182159 -0.18976593
##         V202        V203        V204        V205        V206       V207
## 16 0.9256229  0.37484550  0.67757797  0.61710930  0.73031235  0.3635922
## 17 0.7967205 -0.08593941  0.05968666  0.24624825  0.32953262 -0.1304684
## 18 0.7534370 -0.25172234 -0.23921776  0.06679535 -0.01148605 -0.4939041
## 19 0.6614075 -0.29195595 -0.38718796 -0.09132957 -0.30593872 -0.6345291
## 20 0.5396099 -0.11484337 -0.43234444 -0.18851852 -0.26968765 -0.6574974
## 21 0.4102344  0.05648231 -0.55289078 -0.37616920 -0.16718864 -0.5597649
##          V208       V209        V210       V211       V212        V213
## 16  0.3026562  0.7621098  0.22335625 0.12828064 0.67695236  0.49046707
## 17 -0.4175777  0.2878094  0.03155899 0.04874802 0.57343674  0.14898491
## 18 -0.8276558 -0.0701580 -0.08461189 0.03085899 0.35367203 -0.01460838
## 19 -0.8168735 -0.1740627 -0.01273155 0.22218895 0.28437424 -0.01195717
## 20 -0.8103924 -0.3301563  0.02929497 0.25163841 0.17507553 -0.03765678
## 21 -0.6874218 -0.2782841 -0.07547188 0.06828117 0.09671974 -0.04101562
##          V214       V215       V216       V217        V218       V219
## 16  0.2331219  0.1260948  0.5738277  0.1471119 -0.01969147 -0.3151588
## 17 -0.1932812 -0.5146885 -0.0465622 -0.3875027 -0.53296852 -0.7954674
## 18 -0.3840637 -0.9417210 -0.4767170 -0.7269554 -0.49648476 -0.8614063
## 19 -0.1685925 -0.6757069 -0.2846870 -0.5325794 -0.13094139 -0.6257820
## 20 -0.1003914 -0.3348427 -0.1773453 -0.4710197 -0.06593704 -0.5387478
## 21 -0.1497650 -0.3322658 -0.4166412 -0.5811691 -0.15468788 -0.5660152
##           V220      V221       V222      V223      V224      V225      V226
## 16  0.11265945 0.7871094 0.18210983 0.6070309 0.3407040 0.5929680 1.3456230
## 17 -0.18507957 0.4715614 0.07031059 0.3962479 0.1459370 0.2514858 1.1279716
## 18 -0.29140663 0.3098431 0.12163925 0.4771080 0.1499233 0.1228924 0.9871864
## 19 -0.14812088 0.4196873 0.26851654 0.6059380 0.3931236 0.2492924 1.0926571
## 20 -0.05539131 0.5448437 0.30429459 0.7591400 0.8350773 0.6285915 1.2883587
## 21 -0.22867203 0.4192162 0.22702980 0.7982826 1.0704689 0.8977356 1.4464855
##         V227      V228        V229       V230       V231        V232      V233
## 16 0.3835945 0.4400787  0.23710823 -0.4234409 -0.6201553  0.49265671 0.7646084
## 17 0.6715622 0.4509373 -0.26375008 -0.9267159 -0.9779682  0.07742119 0.4590607
## 18 0.9320297 0.3357830 -0.43195343 -0.8952370 -0.8676548 -0.09015846 0.2360935
## 19 1.1880436 0.3053131 -0.28133011 -0.6159382 -0.6045322  0.02437592 0.3471889
## 20 1.4201565 0.3601570 -0.25977135 -0.5934372 -0.7362499 -0.07538986 0.3373451
## 21 1.5002327 0.6358604 -0.03116989 -0.4821873 -0.7747669 -0.07367325 0.3992157
##           V234        V235        V236         V237       V238        V239
## 16  0.36210823 -0.10046959 -0.22554779 -0.773281100 -0.6468754 -0.34015465
## 17 -0.03343773 -0.38250160 -0.48406219 -1.066015200 -1.0020313 -0.64093590
## 18 -0.10336113 -0.26414108 -0.39007568 -0.840858460 -0.8053150 -0.43421936
## 19  0.11226654  0.04843903 -0.01687622 -0.300706860 -0.3085938 -0.03320503
## 20  0.11304474  0.13039017  0.27382660  0.001092911 -0.1491413  0.03390503
## 21  0.13828087  0.21203232  0.44046974  0.072734833 -0.1610146 -0.08726502
##          V240       V241        V242         V243        V244       V245
## 16 -0.4949207 -0.3678913  0.09031105 -0.381408690  0.08140755 0.21586037
## 17 -0.9828110 -1.0800018 -0.52671814 -0.732967380 -0.08757973 0.14530945
## 18 -1.0229645 -1.3569527 -0.57023621 -0.492654800 -0.01640701 0.04109383
## 19 -0.8709374 -1.1675797 -0.16843987 -0.007030487  0.13937759 0.02593803
## 20 -0.9298458 -1.0422707  0.01156044  0.127500530  0.16585922 0.16359520
## 21 -1.0541401 -0.9186726  0.12406158  0.223983760  0.25757790 0.41796684
##           V246       V247       V248       V249       V250       V251      V252
## 16  0.34710693 0.64578056 -0.1067963 -0.1845303  0.2131214 -0.3064060 0.5575790
## 17  0.01656151 0.21499825 -0.4565620 -0.6985149 -0.4407787 -1.0421886 0.1659374
## 18 -0.32461166 0.03210831 -0.5163288 -0.7821102 -0.6190624 -1.1242218 0.2407837
## 19 -0.33648300 0.09718704 -0.3043766 -0.5294571 -0.5198422 -0.8394566 0.4990635
## 20 -0.15070343 0.13164139 -0.3011723 -0.3976574 -0.4916420 -0.5148430 0.6739044
## 21 -0.13921928 0.08828354 -0.2282810 -0.3597641 -0.5772648 -0.5585155 0.5121098
##         V253        V254       V255       V256        V257         V258
## 16 0.6883602  0.52280807  0.1173420 -0.3673420  0.04461098  0.555856700
## 17 0.3287487  0.31078148 -0.3504696 -0.6413288  0.05531120  0.230312350
## 18 0.2892971  0.30726242 -0.5126534 -0.7239075 -0.05515671  0.002887726
## 19 0.3824196  0.26655960 -0.6420307 -0.6206226 -0.09406090 -0.083982468
## 20 0.3589802  0.09031296 -0.9037476 -0.6853905 -0.25265503 -0.118204120
## 21 0.2575779 -0.02968788 -1.1035175 -0.8199196 -0.34328461 -0.152969360
##         V259         V260        V261        V262       V263       V264
## 16 0.3945312  0.385702130  0.74921799  0.65187263 -0.1651535 -0.5974217
## 17 0.1187477 -0.029062271  0.26773453  0.37172127 -0.1821861 -0.7328110
## 18 0.1196079 -0.096326828  0.13914108  0.29468536 -0.3179703 -0.9342156
## 19 0.1959381  0.005624771  0.04554367  0.17140770 -0.4232063 -0.9621849
## 20 0.2053909 -0.008672714 -0.13515854 -0.00164032 -0.3598423 -0.9435921
## 21 0.1345310 -0.067031860 -0.07351494  0.04273605 -0.3147659 -1.0328903
##          V265         V266       V267        V268        V269         V270
## 16 -0.2053909  0.006557465 -0.4001579 -0.07859039  0.63336182  0.638355260
## 17 -0.3000011 -0.175468440 -0.7092209 -0.43757820  0.20406151 -0.005937576
## 18 -0.4544544 -0.065237045 -0.6789055 -0.42765617 -0.06265831 -0.500860210
## 19 -0.4400787  0.131561280 -0.5432816 -0.34187317 -0.23531342 -0.660232540
## 20 -0.3622704  0.125312810 -0.7737465 -0.52038956 -0.39390564 -0.603204730
## 21 -0.2936726  0.071561813 -1.0147648 -0.67117310 -0.42453384 -0.462970730
##          V271      V272       V273       V274       V275        V276       V277
## 16 0.84953117 0.7444515  1.0642185  0.7081222  0.2785950 -0.02117157 0.24335861
## 17 0.60374832 0.5446892  0.6114845  0.1242218 -0.4334373 -0.40781212 0.29999924
## 18 0.22085762 0.2886734  0.1203918 -0.2240620 -0.7629700 -0.70046616 0.47054482
## 19 0.04843903 0.3018723 -0.2057056 -0.5110931 -0.8394566 -0.84718704 0.55867004
## 20 0.03413963 0.3513260 -0.3414059 -0.8741417 -0.7823429 -0.87984467 0.35023117
## 21 0.04578209 0.2629700 -0.3422661 -0.8860149 -0.5047665 -0.80289078 0.06883049
##         V278        V279         V280         V281       V282      V283
## 16 0.1703110 -0.76391029 -0.133592610  0.282110210 0.59835815 1.1932812
## 17 0.5270328 -0.33671761 -0.086328506  0.086561203 0.25656319 0.8962479
## 18 0.8060112  0.06734657  0.007341385  0.069843292 0.26788902 1.0121078
## 19 0.7203102  0.02672005 -0.040622711  0.004686356 0.26101685 1.0909367
## 20 0.3765640 -0.29874802 -0.186639790 -0.167654040 0.15929413 1.1091404
## 21 0.1015606 -0.49976540 -0.241174700 -0.243284230 0.02453041 1.0557804
##        V284      V285       V286        V287        V288      V289        V290
## 16 1.634453 1.1504688 0.60812187  0.07484436 -0.01992035 0.5021095 -0.02719116
## 17 1.203436 0.5914841 0.08797073 -0.57593536 -0.61906242 0.2437477 -0.38671684
## 18 1.213673 0.6216392 0.12218857 -0.71796989 -0.85421753 0.2405472 -0.51648712
## 19 1.271874 0.6980419 0.38390732 -0.40945625 -0.71343613 0.5386696 -0.34093857
## 20 1.145077 0.5698433 0.44335938 -0.06484222 -0.51359558 0.7352295 -0.25093651
## 21 1.000469 0.4539852 0.40023613  0.03898239 -0.39163971 0.7938290 -0.18718910
##          V291       V292        V293       V294         V295       V296
## 16 -0.4151573 -0.4085922  0.10335922 0.79585648  0.913280490 0.58570290
## 17 -0.6767178 -0.7700787 -0.34844017 0.48280907  0.487499240 0.18218803
## 18 -0.6276531 -0.9126549 -0.44265747 0.25413895  0.223356250 0.06117058
## 19 -0.4332790 -0.8231239 -0.18906403 0.14476585  0.028438568 0.13812065
## 20 -0.3974953 -0.6878910 -0.05265427 0.17429543  0.006639481 0.19507599
## 21 -0.3435135 -0.6411705 -0.18328476 0.06202888 -0.096719742 0.13296890
##           V297        V298        V299        V300        V301        V302
## 16  0.23546791 -0.09062576  0.38109589  0.18008041  0.21461105 -0.06093979
## 17 -0.09226608  0.02172089  0.21531296 -0.11156464  0.08499908 -0.05421829
## 18 -0.11460876  0.03718567  0.08827972 -0.13046646 -0.08820343  0.00601387
## 19 -0.03070831 -0.10734367 -0.05195618 -0.08093643 -0.04382896  0.07156181
## 20 -0.05515671 -0.17663956  0.10515785 -0.05109405  0.30648041  0.20906258
## 21 -0.10726547 -0.19726753  0.10523415 -0.25164032  0.42882919  0.30406380
##           V303        V304      V305       V306       V307      V308      V309
## 16  0.29359055  0.29015541 0.7721100 0.75460815 0.70828056 0.8419533 0.4429684
## 17  0.05327988  0.02116966 0.5815601 0.66781044 0.70374870 0.7959385 0.7502346
## 18 -0.33140373 -0.27640724 0.3110924 0.45788956 0.54585838 0.7011719 1.0791397
## 19 -0.54827881 -0.31187248 0.1621876 0.23226547 0.25843811 0.5168724 0.8530445
## 20 -0.41749954 -0.19038963 0.2785950 0.20304489 0.10039139 0.5563278 0.6723423
## 21 -0.29351616 -0.10367203 0.2454662 0.03703117 0.09703064 0.3904686 0.3277340
##          V310      V311      V312     V313      V314        V315       V316
## 16 -0.1281261 0.1448460 0.5113296 1.208361 0.3628101 -0.19390678 0.65640640
## 17  0.2329693 0.4465618 0.7159386 1.297499 0.4120312 -0.48546982 0.19867325
## 18  0.6484375 1.0332794 1.1445332 1.378048 0.7472629 -0.40265274 0.09484482
## 19  0.4576569 1.0342960 1.1665630 1.243673 0.8503094 -0.17453003 0.38437843
## 20  0.3858585 1.1189060 1.2851543 1.128979 1.1790619  0.08875275 0.61460876
## 21  0.2927342 1.0514832 1.3158588 1.017578 1.2753105 -0.02476311 0.25757790
##          V317      V318      V319      V320       V321       V322        V323
## 16 0.59710884 1.1371059 0.6007824 0.7444515 -0.1945305 0.01437187  0.13234520
## 17 0.35156059 1.0865612 0.7849960 0.7121887  0.1389847 0.21921730  0.05531120
## 18 0.06984329 0.8091392 0.8708592 0.7199230  0.5978908 0.38843536 -0.19922256
## 19 0.09718895 0.4385166 0.5584393 0.5481224  0.6417923 0.31515694 -0.35320473
## 20 0.28609276 0.3030453 0.3666420 0.4713268  0.6698418 0.49085999 -0.07734299
## 21 0.11046600 0.1432781 0.2645302 0.3579712  0.5952339 0.47398567 -0.16976547
##           V324       V325       V326       V327       V328
## 16  0.43882942 -0.3291397  0.7328091  0.8723450  0.4501572
## 17  0.23093796 -0.7037506  0.2107830  0.5207806  0.2324200
## 18 -0.06046677 -1.2319546 -0.6264858 -0.1476555 -0.1301556
## 19 -0.04093361 -1.1313286 -0.6496925 -0.2182808 -0.0706253
## 20  0.04640579 -0.7122688 -0.5821876 -0.3124962 -0.1203918
## 21 -0.22789192 -0.8699207 -0.7609367 -0.5010166 -0.3286724
t(SSTdata) %>% dim()
## [1]  328 2261
# construct the EOFs 
Z <- t(SSTdata) # data matrix, rows are times, columns are spatial sampling locations.
spat_mean <- apply(SSTdata, 1, mean) # average across all time.
nT <- ncol(SSTdata) # no of time points
# outer command makes the spatial mean into a matrix
Zspat_detrend <- Z - outer(rep(1, nT), spat_mean) # detrend data (center to spatial mean)
Zt <- 1/sqrt(nT-1)*Zspat_detrend # normalize
E <- svd(Zt) # SVD
# number of EOFs to estimate in the data
n <- 10

# Projected onto the 
TS <- Zspat_detrend %*% E$v[, 1:n]
summary(colMeans(TS))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2.091e-16 -8.031e-17  3.201e-16  3.113e-16  6.916e-16  8.885e-16

This verifies that the time series have a zero mean.

tau <- 6
nT <- nrow(TS)
TStplustau <- TS[-(1:tau), ] # TS with first tau time pts removed 
TSt <- TS[-((nT-5):nT), ] # TS with last tau time pts removed
Cov0 <- crossprod(TS)/nT # cross product
Covtau <- crossprod(TStplustau, TSt) / (nT - tau)
C0inv <- solve(Cov0)
Mest <- Covtau %*% C0inv # estimated M, the dynamics
Ceta <- Cov0 - Covtau %*% C0inv %*% t(Covtau) # Covariance matrix
image(Mest)

image(Ceta)

The forecast is given by \(\hat\mu + \PhiM^2\alpha_t\).

SSTlonlat$pred <- NA
alpha_forecast <- Mest %*% TS[328, ] # propigate alpha first,
idx <- which(SSTlonlat$mask == 0)
SSTlonlat$curr[idx] <- as.numeric(E$v[, 1:n] %*% TS[328, ] + spat_mean)
SSTlonlat$pred[idx] <- as.numeric(E$v[, 1:n] %*% alpha_forecast + spat_mean) # add to spatial mean for original scale.

# for plotting
SSTlonlat$obs1[idx]  <- SSTdata[, 328]
SSTlonlat$obs2[idx]  <- SST_Oct97


C <- Mest %*% Cov0 %*% t(Mest) + Ceta # 6 month ahead variances

SSTlonlat$predse[idx] <- sqrt(diag(E$v[, 1:n] %*% C %*% t(E$v[, 1:n])))

SSTlonlat %>% ggplot(aes(lon, lat, fill = pred)) +
  geom_tile()

See the appendix for proper visualization.

DSTM_EM is provided with the package STRbook, and runs an EM algorithm that does ML estimation in state-space models.

DSTM_Results <- DSTM_EM(Z = SSTdata,
                        Cov0 = Cov0,
                        muinit = matrix(0, n, 1),
                        M = Mest,
                        Ceta = Ceta,
                        sigma2_eps = 0.1,
                        H = H <- E$v[, 1:n],
                        itermax = 10,
                        tol = 1)

The object DSTM_Results has the setimated parameters, and complete-data negative log-likelihood. Estimates of \(\alpha_{1,t}\), \(\alpha_{2,t}\) and \(\alpha_{3,t}\)

par(mfrow = c(1,3))
for(i in 1:3) {
  plot(DSTM_Results$alpha_smooth[i, ], type = "l",
       xlab = "t", ylab = bquote(alpha[.(i)]))
  lines(TS[, i], lty = "dashed", col = 'red')
}

image(as(DSTM_Results$Mest, "Matrix"))

image(as(DSTM_Results$Mest %^% 6, "Matrix"))

image(as(Mest, "Matrix"))

alpha <- DSTM_Results$alpha_smooth[, nT]
P <- DSTM_Results$Cov0
for(t in 1:6) {
  alpha <- DSTM_Results$Mest %*% alpha
  P <- DSTM_Results$Mest %*% P %*% t(DSTM_Results$Mest) + DSTM_Results$Ceta
}